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Some Kondo insulators are expected to possess topologically protected surface states with linear 
Dirac spectrum, the topological Kondo insulators. Because the bulk states of these systems typically 
have heavy effective electron masses, the surface states may exhibit extraordinarily small Fermi 
velocities that could force the effective fine structure constant of the surface states into the strong 
coupling regime. Using a tight-binding model we study the many-body instabilities of these systems 
and identify regions of parameter space in which the system exhibits spin density wave, and charge 
density wave order. 


I. INTRODUCTION 

In many metals and semiconductors the behavior of 
the low energy electronic states can be understood in 
terms of free quasiparticles with quadratic energy dis¬ 
persion in momentum p 2 /2 m*, where m* is a system- 
dependent effective massi. However, there are a number 
of materials whose low energy electronic states are bet¬ 
ter described as massless Dirac fermions, including the 
superfluid phase of 3 He, high-temperature d -wave super¬ 
conductors, graphene, and the surface states of topolog¬ 
ical insulators^ —. In these Dirac materials the kinetic 
energy is proportional to the momentum vp, just like 
massless relativistic particles but with a speed v that 
depends on the details of the system. For example, in 
graphene v « 10 6 m/s« c/300. 

The fact that quasiparticles obey the Dirac equation 
instead of the Schrodinger equation can affect a variety of 
electronic properties, for example, the integer quantum 
Hall effect and localization^. Another important way the 
Dirac nature of the quasiparticles manifests itself is in 
the effect of interactions. If the quasiparticles of a sys¬ 
tem obey the Schrodinger equation, then the ratio of the 
average interparticlc Coulomb energy to the average ki¬ 
netic energy, r s = Ec/Ek , is related to the density by 
r s oc where the constant of proportionality de¬ 

pends on characteristics of the material. In contrast to 
normal metals, for Dirac materials this ratio is a charac¬ 
teristic of the system, independent of the electron density, 
given by a = Ec/Ek = e 2 /(hev). In this expression e 
is the charge of the electron, e is the material’s dielec¬ 
tric constant, h is the reduced Planck constant, and v is 
the speed of the Dirac particles. Much work has gone 
into the study of the phase diagram of graphene with re¬ 
spect to this parameter o£ and the results indicate 
that there is a critical value a c such that if a < a c the 
spectrum remains gapless and if a > a c the system flows 
toward the strong coupling regime and is likely to develop 
a gap&. Thus far, perturbative and numerical results sug¬ 
gest the critical value is a c « 1& & 12 i 13 , while experiments 


involving suspended graphene, for which a ~ 2.2, seem 
to indicate a gapless state to within 0.1 meV of the Dirac 
point^. Therefore, the ground state of Dirac materials in 
the strong coupling regime is not currently understood. 
For this reason we propose studying a class of materials 
with much smaller Fermi velocity than that of graphene 
since this class of materials is likely to possess a a c 
and would be a better candidate for experiments probing 
the strong coupling regime in Dirac materials. 

The surface of a three-dimensional (3D) topological in¬ 
sulator (TI) hosts two-dimensional (2D) Dirac quasipar¬ 
ticles similar to those found in graphene. Examples of 
experimentally verified 3D TIs include Bi 2 Se 3 , Bi 2 Te 3 , 
and Sb 2 Te 3 , all of which have Fermi velocities roughly 
half of that in graphen o 15 ! 16 . However, there is another 
class of topological insulators, the topological Kondo in¬ 
sulators (TKI), in which the bulk states are formed by 
renormalized /-electron levels which hybridize with con¬ 
duction electrons to form a millivolt-scale gap in the bulk 
spectrum~~— . The small gap in these materials com¬ 
bined with the large bulk effective mass imply that the 
surface Fermi velocity could be quite small. Some ma¬ 
terials theoretically predicted to fall into this category 
include SmBg^i, YbBi^, and PuBg2£. Furthermore, 
there is a growing body of experimental evidence demon¬ 
strating that SmBg does in fact host metallic surface 
states^—. 

Previous work has explored the possibility of broken 
symmetry states on the surface of a TKI employing a 
continuum model^i. In this paper we present a tight- 
binding model to study the surface states of a TKI and 
proceed to investigate the possible ordered ground states 
for these systems within a mean field theory. From this 
analysis, we find regions of parameter space for the model 
that admit spin density wave and charge density wave 
solutions. For the case of strictly repulsive interactions 
we find that these ordered solutions lie within the region 
of parameter space corresponding to the strong coupling 
regime of Dirac materials (a > a c ps 1). 
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II. THEORETICAL MODEL AND METHODS 


To model the 2D surface states of a TKI, we consider 
a Hamiltonian defined on a square lattice: 

H ° = ' C&J x <T ) 

^ a TV V '“Z 3 

a.P.o- (ij) /-g 

a,yS.cr 


where a and /? are orbital indices, u is a spin index, Rij 
is the unit vector pointing from lattice site j to lattice 
site i, and the matrix T is defined as 

(4 T ; i = j 

Ttj = < —Y ; i,j nearest neighbors . (2) 

I 0 ; otherwise 


The term proportional to A leads to the formation of 
four separate Dirac points in the Brillouin zone. The 
term proportional to T acts as a momentum-dependent 
mass term which gaps out all of the Dirac points except 
the one at k = 0 allowing the model to represent the 
surface states of a strong TI^. The energy eigenvalues 
associated with this Hamiltonian in k space are given by: 


Et = ±4 r 


ak T 


Sill 


sin' 


2 Otky 


\ 


. . A \ 2 sin 2 ak x + sin 2 ak v ^ 


sm 


■ Sill' 


2 ak v 


Expanding this dispersion for small k along the k x direc¬ 
tion we find: 


ft: ± ( aAk + 


(3T 2 - A 2 )c 

6 A 


( 4 ) 


Thus, we can see that to first order in k the dispersion 
matches the Dirac dispersion with Fermi velocity given 
by aA/h. In Fig. [T]we plot the full dispersion from Eq. 
([3j for different ranges of k to demonstrate the Dirac dis¬ 
persion for a few different values of the Fermi velocity. It 
shows that near the Dirac point the parameter A con¬ 
trols the Fermi velocity; however, for 4<Fwe can see 
that the cubic term in Eq. £]} begins to dominate and 
the dispersion away from the Dirac point becomes no¬ 
ticeably less linear. Since we are most interested in the 
regime in which the model best describes a Dirac ma¬ 
terial, in this paper we focus on the case in which the 
chemical potential is close to the Dirac point. 

In terms of the model parameters the bandwidth is 
given by: 


w 


16 r 

4 A 2 

y/2A 2 + 16T 2 


; A < 4T 
; A > 4T 



FIG. 1. (Color online) Plots of the band structure given by 
Eq. m along the diagonal of the square Brillouin zone using 
four different values of the parameter A [1/4 (solid, black), 
1/8 (dashed, red), 1/16 (dashed-dot, blue), and 1/32 (dot¬ 
ted, green)], in two different momentum ranges [(a) from 
k = (—l/5a, — l/5a) to k = (l/5a, l/5a) and (b) from 
k = (— n/a, — Tt/a) to k = (n / a,n / a)\. All energies are in 
units of the bandwidth. 


Note that for A < AT the bandwidth is a constant set by 
the model parameter T. In the analysis that follows we 
restrict the range of A to A < 4T and present all energies 
in units of the bandwidth w = 16T. We also present all 
distances in units of the lattice constant a. 

To account for interactions we consider the full Hamil¬ 
tonian: H = Hq+Hi, where Hi takes on the exact form: 




0 -l r «— r dA 


a, a' i^j 


V / irT r r 7 F+^ 




- u 55 

i 

( 5 ) 

where Vg controls the strength of the long-range Coulomb 
interaction between / electrons and U is introduced as 
an on-site interaction between / electrons. 


In our calculations, we replace the exact interaction 
term Hi with the mean field Hamiltonian: 


o— |r;-r,-|/A 




+ (a iV’lj.tV’j,/,; + t) + E 0- 

( 6 ) 

To capture the on-site Coulomb repulsion between / elec¬ 
trons we can set U = —Vo/d. If we wish to include an 
attractive on-site interaction we set U > 0. In the ab¬ 
sence of a large electron-phonon coupling this attractive 
interaction could be engineered by the adsorption of a 
finite density of nonmagnetic molecules as discussed pre¬ 
viously in the context of topological insulators^. We 
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may choose to write this in a more compact notation as: 

Hj IF = g a ,1pj t p,a' 

a,(3,cr,(j' i,j 

+ Y t) + E 0 

i 

where 

y ° J\ ri _ r . 3 \2 +d 2{m,f,*)8afi6 a f ; i^j 
U (jli^ f^a^Sczfi &af (1 ^erer') ) ^ — J 


where /(e) = et y fcj ) r+1 is the Fermi-Dirac distribution 
function at temperature T and ks is the Boltzmann con¬ 
stant. Given an initial set of model parameters and 
a temperature, Eqs. © and © allow us to solve for 
the density profile and superconducting order parameter 
A self-consistently. In the next section we discuss our 
progress toward solving these equations. 

In some cases multiple solutions for the same model 
parameters may be found. In this case it is useful to 
compare the free energy associated with each of the so¬ 
lutions, given by F = ksTlnZ, where Z is the partition 
function. The true ground state of the system will be 
given by the solution with the lowest free energy. 


A* — /,tV’i,/,!)• 

Equipped with this mean-field Hamiltonian we perform 
a Bogoliubov transformation: 

rfitot ,f = — lni V i,ot,n,\ 

n 

= y ( r ynl u i,a,n,l + lf n ^ v i,a,n,i^ 
n 

where 7 ^ {"/na) creates (annihilates) an eigenstate of 
the mean-field Hamiltonian H. It can be shown that the 
coefficients u and v satisfy the following equations: 

Cn7^ J c*,n,t' = y ( E ia.'\,j 0 '\ r ^ l j, 0 ,n,'\ A 
3,0 

Cn,t v i,u,n,l = ~ 'Y, E ia\.,j0X v 3,0,n,i + ^i u i,f,n,t 

(7) 

4.^2,a,n,4. — / J T A iVi,f^ n ^ 

3,0 

e n,l v i,a,n,t = — y ] E iat,j0t V j,0,n,t + 

3,0 

where Hiaa,j E ioa.jfiu' Wi&a.jj3a ! and e n cr are 

eigenvalues of H. 

Given the solutions to these equations we can write the 
mean fields as 
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n 
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III. NUMERICAL RESULTS AND DISCUSSION 

While it is straightforward to numerically solve Eqs. 
© and © for a finite system, we can make the compu¬ 
tation more efficient by using the supercell technique^ 
as described in the appendix. For a system with a 10 x 10 
real space unit cell and an 8 x 8 supercell we solved for 
self-consistent solutions to Eqs. © and ©. For the 
following results we focused on the case of nearest neigh¬ 
bor Coulomb interactions only and the zero temperature 
limit. In Eq © we used a screening length of A = 1 and a 
lattice cutoff of d = 1. We considered two limiting cases: 
the case of a repulsive on-site interaction (U = —Vo), and 
the case of an attractive on-site interaction that scales 
with the Coulomb interaction {U = Vo). 

Starting from initial seeds that possessed antiferromag¬ 
netic, ferromagnetic, checkerboard and stripe charge den¬ 
sity wave (CDW) order in addition to random seeds we 
found self-consistent solutions for Eqs. © and © us¬ 
ing a convergence criterion of 10” 3 . Some of the self- 
consistent solutions that emerged from the different seeds 
for the same model parameters differed from each other. 
In these cases the one with the lowest free energy was 
taken to be the solution. In Fig. [2] we show the regions 
of parameter space for which we found solutions in the 
case of repulsive on-site interactions while in Fig. [4] we 
show the regions of parameter space for which we found 
solutions in the case with attractive on-site interactions. 

First, we consider the case of on-site repulsion (U = 
—Vo), Fig. [21 Note that the general trend is consistent 
with our expectations for Dirac materials. In the region 
of strong coupling, a = Vq/A > a c , we find Coulomb- 
driven ordered states, while in the weak coupling region, 
Vo/A < a c , a paramagnetic (PM) normal metallic state 
exists. These results are consistent with the established 
value of a c « 1. However, it appears that there is a 
critical value of the coupling, V c ~ w/3, for this model 
below which the solution is trivial. This is in contrast 
to the case of a Dirac continuum model in which the 
only parameter governing the Coulomb interaction is a. 
This difference can be attributed to the fact that for very 
small values of A the band structure appears less linear 
and eventually the cubic term becomes more important, 
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FIG. 2. (Color online) Plot of the phases for the self-consistent 
solutions found in different regions of the A,Vo plane. Note 
that the region below the line A = Vo appears to favor the 
formation of nontrivial order, consistent with a c ~ 1. The 
region enclosed by the red dashed line favors the formation of 
spin density wave order while in the region enclosed by the 
black solid line we find both spin density wave and charge 
density wave solutions. Outside of these regions the solution 
is paramagnetic (PM). 


as we can see from Eq. ([4]). It is reasonable to expect that 
real materials which host slow Dirac states will typically 
have similar behavior since the bands for these materials 
are expected to develop nonzero curvature away from the 
Dirac poin t 18 d 9 i 24 i 29 . 



FIG. 3. (Color online) Plot of the density modulations over a 
10x10 real space unit cell, as observed in the SDW and CDW 
regions shown in Figs. [2] and [4] 

Taking a closer look at Fig. [2]we can see that there are 
three distinct regions of the Vo, A plane: a region favoring 
spin density wave (SDW) order, a region in which SDW 
and CDW coexist, and a region in which the solution was 
PM. Both SDW and CDW modulations were associated 


with ( 7 r, 7r) wave vectors as shown in the sample plot in 
Fig. [3j Intermediate states were also observed but these 
appear to be higher energy excitations. In the SDW re¬ 
gion the boundary for the phase along the Vo axis, at 
approximately one third of the bandwidth, defines the 
critical coupling, V c . We find that above another critical 
value of Vo CDW order begins to coexist with the SDW. 
In the coexistence region for some model parameters we 
were able to find solutions with exclusively CDW order 
but we lacked the resolution to see if these solutions in¬ 
dicated the existence of an additional region of the plane 
in which CDW is truly favored over SDW order, further 
calculations will be needed to answer this question. 

Next we turn our attention to the case in which we 
include an on-site attraction (U = Vo), as shown in Fig. 
U In this case we find two regions: a region with CDW 
and a PM region. Again, these density modulations were 
associated with a wave vector of (7r,7r). The region that 
favors CDW begins at Vo ~ w /3 and covers the rest of 
the plane. It is interesting to note that the CDW or¬ 
der appears for Vo > w /3 which is the same as V c for 
the case with repulsive on-site interactions. It should be 
noted that some of the self-consistent solutions we found 
near the transition region Vo ~ w/3 seemed to possess 
a small superconducting order parameter; however, this 
order parameter was usually just below the convergence 
criterion (even when the convergence criterion was low¬ 
ered to 10- 7 ). We attribute the absence of a supercon¬ 
ducting region to the fact that we restricted ourselves 
to the case of half filling in which there was no density 
of states to allow for superconducting pairing. A more 
detailed study of the region near Vo « w/3 may be in¬ 
teresting for future work studying this model away from 
half filling. 

Note in Fig. 0]the absence of any regions with magnetic 
order, in contrast to Fig. [2] in which both AFM and FM 
order were found. This can be accounted for by a heuris¬ 
tic argument based on Eq. ©• Notice that the spin- 
dependent terms in the mean field Hamiltonian are given 

by ~ U E i > thus 

the expectation value of the contribution to the total 
energy will be -217 Ei( n i,/,t)( n b/,4.)- For U > 0 we 
can see that the energy can be minimized if the sum 
Ei( n i,/,t)( n b/4) t a kes on its maximum possible value. 
Each term of this sum has a maximum value when 
( n iJ, t) = = 1/2- Therefore the minimum energy 

can be expected to be achieved in a state with no mag¬ 
netic order. However, for U < 0 the minimum energy is 
achieved for a minimum value of Ei( n i,/,t)( n * /,4.)j which 
can allow the system to minimize its energy through an 
on-site spin polarization. 


IV. CONCLUSIONS 

In conclusion, we presented a model for studying the 
surface states of a class of topological Kondo insulators 
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FIG. 4. Plot of the phases for the self-consistent solutions 
found in different regions of the A,Vo plane for attractive on¬ 
site interaction. In the region enclosed to the right of the 
solid black line the self-consistent solutions possessed charge 
density wave order; outside of this region the solution was 
paramagnetic (PM). 
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Appendix A: Supercell Technique 

The system of equations given by Eqs. 0 and (0 can 
be solved for finite systems by simple matrix diagonal- 
ization. However, the matrix that must be diagonalized 
is 8 N x 8 N, where N is the number of lattice sites and 
8 = 2(spins) x 2(orbitals) x 2(electron-hole). We can 
see that for 6400 sites this would involve diagonalizing 
a 51200 x 51200 matrix which is not terribly practical. 
Using the supercell technique we can decrease the size 
of the matrix that needs to be diagonalized significantly. 
In the framework of the supercell technique we recognize 
that, due to the periodicity of the system, the solutions 
u ri , a ,n,c t and u riiQ>rii0 ., are Bloch waves. To account for 
this we write 

j.OL.n.cy — ^ 1 ^k,ri,a,n,cr , . 

Ai) 

i,a,n,(T ^ ^k,r i,a,n,cn 

where k is the crystal momentum. After this transfor¬ 
mation Eqs. 0 and @ become: 


and explored the dependence of the band structure on 
the model parameters, identifying the parameters which 
determine the Fermi velocity at the Dirac point. We 
then added interactions to this model, accounting for 
both Coulomb interactions as well as the possibility of 
an on-site attractive interaction. Using mean-field the¬ 
ory, at zero temperature, we found self-consistent solu¬ 
tions for different model parameters, investigating the re¬ 
lationship between the Fermi velocity at the Dirac point, 
the strength of the interactions, and the nature of the 
self-consistent solutions. For the case with on-site repul¬ 
sion we identified three regions of parameter space with 
different Fermi velocity and coupling strength: a region 
which exclusively favored spin density wave order, a re¬ 
gion of coexisting spin density wave and charge density 
wave order, and a paramagnetic normal metallic region. 
We also identified a critical value of the Coulomb in¬ 
teraction strength V c ss w/3 below which the solutions 
were normal metallic. When we considered the case of 
an attractive on-site interaction we found that the so¬ 
lutions possessed charge density wave order above this 
same critical Coulomb interaction strength. 
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where k = nfe where = 1,2,..., M x and 

n y = 1,2,..., M y , M x and M y are the number of unit cells 
in the x and y direction, respectively, M xy = M x M y , fV x 
and IVy are the number of lattice sites per unit cell in the 
x and y direction, respectively, and we define 


Hi, 


ioL<j,jfi<j' ;k 


E y e i k - (r j+Ri- r i) H 


R 


a:cr,(iv 7 -+Rj)/3cr' • 


Now, a system composed of 6400 sites can be studied 
by diagonalizing a 10 x 10 real space system using an 
8x8 supercell. This means we only need to diagonal¬ 
ize a 800x800 matrix instead of 51200 x 51200. More¬ 
over, this diagonalization is performed for each k inde¬ 
pendently and thus the procedure may be easily paral¬ 
lelized to further improve performance. 
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